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Abstract 

We consider molecular communication networks consisting of transmitters and receivers 
distributed in a fluidic medium. In such networks, a transmitter sends one or more signalling 
molecules, which are diffused over the medium, to the receiver to realise the communication. 
In order to be able to engineer synthetic molecular communication networks, mathematical 
models for these networks are required. This paper proposes a new stochastic model for 
molecular communication networks called reaction-diffusion master equation with exogenous 
input (RDMEX). The key idea behind RDMEX is to model the transmitters as time series 
of signalling molecule counts, while diffusion in the medium and chemical reactions at the 
receivers are modelled as Markov processes using master equation. An advantage of RDMEX 
is that it can readily be used to model molecular communication networks with multiple 
transmitters and receivers. For the case where the reaction kinetics at the receivers is linear, 
we show how RDMEX can be used to determine the mean and covariance of the receiver 
output signals, and derive closed- form expressions for the mean receiver output signal of the 
RDMEX model. These closed-form expressions reveal that the output signal of a receiver 
can be affected by the presence of other receivers. Numerical examples are provided to 
demonstrate the properties of the model. 

Keywords: Molecular communication networks, nano communication networks, synthetic molec- 
ular communication networks, master equations, stochastic models, synthetic biology 

1 Introduction 

We consider molecular communication networks consisting of transmitters and receivers dis- 
tributed in a fluidic medium. In such networks, a transmitter sends one or more signalling 
molecules, which are diffused over the medium, to the receiver to realise the communication. 
The study of molecular communication has its origin in biology and biophysics. Molecular com- 
munication is a vital mechanism in multi-cellular organisms. The human body, which has an 
estimated 10^^ cells, uses molecular communication to keep the body in a healthy state. In fact, 
cells in the human body constantly communicate with other cells using molecular communication. 

There are a couple of reasons why synthetic molecular communication networks, which are 
inspired by molecular communication in living organisms, should be studied [H EJ [3]. Firstly, 
synthetic molecular communication networks can be combined with nano-sensors and molecular 
computing [4 to form nano-sensor networks [5] for health monitoring, medical diagnosis and 
cancer therapy. Secondly, the study of synthetic molecular communication networks can be used 
to enhance our understanding of their biological counterparts. 

In order to be able to engineer synthetic molecular communication networks, we need math- 
ematical models which can be used to predict the performance of these networks. For example, 
if a transmitter in a molecular communication network emits a number of signalling molecules 
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to communicate with a receiver, we would like to be able to determine the receiver output sig- 
nal in order to determine the probability of correct reception at the receiver. Such evaluations 
can be realised if a mathematical model is available to determine receiver output signal based 
on the transmitter's input signal. The main contribution of this paper is that we propose a new 
stochastic model for molecular communication networks. Our model is based on reaction-diffusion 
master equation (RDME) [6 which is a well known model in physics and chemistry for modelling 
systems with both diffusion and chemical reactions. In this paper, we propose an extension to 
RDME, which we call reaction-diffusion master equation with exogenous input (RDMEX). The 
key idea behind RDMEX is to model the transmitters as time series of signalling molecule counts, 
while diffusion in the medium and chemical reactions at the receivers are modelled by RDME. 
An advantage of RDMEX is that it can readily be used to model molecular communication net- 
works with many transmitters and many receivers. For the case where the reaction kinetics at 
the receivers is linear, we show how RDMEX can be used to determine the mean and covariance 
of receiver output signal of molecular communication networks. These results allow us to derive 
closed-form expressions showing how the receiver outputs relate to the transmitter signals when 
there are multiple transmitters and receivers. These expressions show that the output of a receiver 
can be influenced by the presence of other receivers in a molecular communication network. They 
also reveal the coupling between diffusion and chemical reactions at the receivers. 

This paper is organised as follows. In Section [2| we present background materials on master 
equations. In Section [3j we present the RDMEX and show how it can be used to model molecular 
communication networks with multiple transmitters and receivers. We also show in this section, 
for the case where the reaction kinetics at receivers is linear, how we can determine the mean 
and covariance of the receiver output signals in molecular communication networks. The rest 
of the paper is focused on determining the mean receiver output signal. We approach this by 
using two methods, which will be discussed in sections [4] and [Sj In Section [4j we determine the 
continuum limit (i.e. infinite spatial resolution) of the RDMEX and show that it results in a 
reaction-diffusion partial differential equation (RDPDE). We derive a closed-form solution to this 
RDPDE and interpret the results. In Section [5) we determine the mean receiver output signal of 
RDMEX with finite spatial resolution and derive another closed-form solution. Numerical results 
are then presented in section [6] to show the accuracy of our solutions. Finally, section [7| describes 
the related work and section Is] gives the conclusions. 



2 Background on master equations 

The aim of this section is to provide the necessary background on master equations. The treatment 
here is brief and includes only the results needed for this article. The reader can refer to the texts 
[HI [7] or tutorial article [8 for a more complete treatment of this subject. This section is divided 
into two parts. We first introduce the general master equation and give a simple example on how 
to use master equation to model a chemical reaction. We then quote some results on mean and 
covariance of the Markov processes modelled by master equations. 



2.1 General master equation 

Consider a continuous-time integer- value vector Markov process Q{t) G Z^, where p is the number 
of vector components, Z is the set of all integers and t is time. When the Markov process Q{t) is 
in state g G Z^, it jumps to the state q + rj (where rj G Z^ with j = 1, 2, ...J and J is the total 
number of possible jumps) at a transition rate of Wj{q). Let P{q^t\qo^to) denote the conditional 
probability that Q{t) = q given that Q{to) = q^. 

We are interested to determine how P{q^t\qo^to) evolves over time. We can do this by using a 
coupled set of ordinary differential equations (ODEs) known as the master equation: 

dP{q,tJqo.to) ^ ^iy,-(g-r,)P(^-r,,t|^o,to)-^iy,(^)Pfet|^o,t^ (1) 
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where one equation of the form ([T| is needed for each vahd state q. Note that the first and second 
terms on the right-hand side of IT]) can be interpreted, respectively, as the rates of entering and 
leaving the state q. In order to simplify notation, we will write P{q^ t) instead of P{q^ t\qo^to) from 
now on. 

A common application of master equation is to model the dynamics of chemical reactions [9]. 
We will give a simple example to illustrate that. 

Example 1 Consider the chemical reaction: 



L + C 

k- 



where the chemical species are L, R and C , and the forward and reverse reaction constants are A:+ 
and k- respectively. This chemical reaction can he described by a Markov process with state space 
Q = [^L^^R^^c]^ where ul is the number of molecules of chemical L etc., and T denotes matrix 
transpose. 

There are two possible types of jumps (i.e. J = 2) in this system. The forward reaction is 
modelled by the jump ri = [—1, —1, 1]^ where the entries of ri reflects the fact that one molecule 
of L and one molecule of R react to form a molecule of C. The rate of jump, which according 
to standard result in chemical kinetics, is Wi{q) = k-^q{l)q{2) = kj^nLUR where is the 
first component of the vector q, etc. Similarly, the reverse reaction is modelled by the jump 
r2 = [1, 1, —1]"^ with W2{q) = k-q{3) = k-Uc- The master equation for this chemical reaction is: 

2 2 

= Y.Wj{q-r,)P{q-rj,t)-J2Wj{q)P{q,t) (2) 

Equations of the type ^ is also known as chemical master equations. 
2.2 Results on mean and covariance 

The master equation ([T]) shows the time evolution of the state probability of the Markov process 
Q{t). However, ([T]) can be difficult to work with because we need one equation for each valid state, 
hence the number of equations can be very large. Therefore, it is easier if we can determine the 
mean and covariance of the Markov process which is defined as follows: 

{Q{t)) = ^9Prob(Q(i)=9) = ^gPfei) (3) 

q q 

m = Y.^q-{Qm{q-{Q{t))fP{q,t) (4) 

q 

where (•) will be used in this paper to denote the mean operator. 

It is possible to use ([T]) to derive how the mean and covariance of the state of the Markov 
process Q{t) evolve over time. The following proposition is taken from ^IQj . 

Proposition 1 For the general master equation we have 



d{Q{t)) 



dt 

In particular, ifWj{q) is a linear function of q. Let Xlj=i ^j^j(^) — then 



dt 
dT,{t) 
dt 



= Am + mA'' + j2^,rjw,{{Qm (7) 

i=i 
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3 Modelling molecular communication network using RD- 
MEX 



We consider a molecular communication network with multiple transmitters and multiple receivers 
in an isotropic fluidic medium, see Figure [l] The communication takes place by the transmitter 
emitting one or more signalling molecules over time. Once the molecules leave the transmitter, 
they diffuse in the medium according to Brownian motion. The receivers are assumed to consist of 
one or more receptors. When a signalling molecule L reaches a receptor they may bind together 
to form a complex C. We will consider the number of complexes at the receiver as the output of 
the receiver. For example, a receiver may infer that a bit has been sent by the transmitter when 
the number of complexes exceeds a threshold. Note that use of this definition of output is also 
used in chemotaxis models in biophysics [TTl [12] and molecular communication network models in 
engineering [13]. 

Based on the above description, we see that a model for molecular communication networks 
must at least capture the diffusion of signalling molecules and the reactions at the receivers. Both 
reactions and diffusion can be modelled as Markov processes, so the master equation ([T]) is a natural 
choice. The remaining issue is how we can model the transmitter. In this paper, we model the 
transmitter by a time sequence which specifies the number of molecules emitted by the transmitter 
at a particular time. This approach is fairly general and can be used to model encoding methods 
that have been considered in the literature, such as molecular coding [1 and concentration coding 
[I4j . We will consider other modelling approaches, e.g. modelling the internal mechanism of the 
transmitters, in future work. 

In section [TH we introduce the RDMEX model by way of an example and then we prove some 
results on mean and covariance of the RDMEX model in section [321 

3.1 The RDMEX model 

In this section, we will introduce the RDMEX model for a molecular communication network 
with 2 transmitters, 2 receivers in a 1-dimensional medium. The reason for that is to simplify the 
presentation. It is fairly simple to generalise the model to multiple transmitters, multiple receivers 
in a 3-dimensional medium, which we will discuss at the end of the section. 

Another simplification is that we will assume, at the receiver, the rate at which the complexes 
are formed is a linear function of the number of local signalling molecules and is independent 
of the number of receptors. (This will be made precise below.) This simplification allows us to 
produce closed form expressions in the continuum limit. Note that, it is straightforward to model 
non-linear reaction rates or incorporate more complex receivers in our model. 

The following is a list of model assumptions, parameters and notation. 

1. We assume that both transmitters use one and the same type of signalling molecule L. 

2. For transmitter 1, we assume that it emits ki^i signalling molecules of L at time ti^i, ki^2 
molecules of L at time ti,2, ^i,6 molecules of L at time ^1,5, where ki^^ (6 = 1, 2, 3, ...) are 
positive integers and ti^^ G R. Similarly, transmitter 2 emits /c2,6 signalling molecules of L at 
time ^2,5 where 6 = 1,2, 3, .... The number of molecules ka^b emitted at time ta^b is assumed 
to be independent of the state of the system at or before ta^b- 

3. The medium is assumed to be a 1-dimensional space of length X. The medium is partitioned 
into TV equal width voxels of width Ax such that N Ax = X. We index the voxels by using 
1,...,A/'. See Figure [2] for an illustration. 

4. The medium is assumed to be isotropic. The rate at which a signalling molecule L diffuses 
from one voxel to a neighbouring voxel is d per molecule per unit time. The rate of diffusion 
from a voxel to a non-neighbouring voxel is zero. Also, the molecule cannot leave the medium 
and we assume the boundary is reflective. The parameter d is related to one-dimensional 
macroscopic diffusion constant by d = 
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5. We assume that each transmitter or receiver occupies only one voxel. Transmitters 1 and 
2 are located respectively in the voxel indexed by Ti and T2. Similarly, we assume that 
receivers 1 and 2 are located in the voxels indexed by Ri and R2. The indices Ti, T2, Ri 
and R2 are integers in the interval [l^N] and are assumed to be distinct. (Note that it is 
simple to modify the model so that a transmitter or a receiver occupies multiple voxels.) 

6. We assume both receivers 1 and 2 use the same type of receptors R and these receptors are 
fixed in space, i.e. do not diffuse. With the simplification mentioned in the introductory 
part of this section, we assume that the reaction between the signalling molecule L and the 
complex C is: 




where /c+ and k- are the macroscopic reaction rate constants. With this assumption, the 
rate of formation of complexes at receiver 1, at any time, is proportional to the number of 
signalling molecules L in the Ri-th voxel. The situation at receiver 2 is similar. 

Remark 1 The reaction kinetics assumed above can be viewed as a linearisation of the 
second order reaction L -\- R ^ C. Similar assumption is also used in fTl\ [T^ to model 
receptor kinetics in chemotaxis. 

7. The state vector q consists of {N + 2) elements where 

q = [ riL,! nL,2 .... til^n nc,i nc,2 ] (8) 

where n^j- represents the number of signalling molecules at the j-th voxel and nc,n represents 
the number of complexes at the u-ih. receiver. 

8. We define two indicator vectors ILti,1t2 ^ Z^^^. The Ti-th element of is 1 and are 
otherwise zero. 1^2 is similarly defined. 

9. In order to write down the diffusion and reaction within the molecular communication net- 
work, we define the following state transition vectors rj and transition rates Wj. The total 
number of possible jumps in this system is J = 2N + 4 where 2N of them model diffusion 
and the rest models reactions at the receivers. We will state these jumps below in four 
categories. Note that all rj G Z^+^ and only the non-zero elements of rj are stated, and q 
is the state vector defined in ([8|. 

(a) The diffusion of L from voxel j to j + 1, where 1 < j < — 1, is modelled by rj and 
Wj{q). Specifically, r^(j) = -1, rj{j + 1) = 1 and Wj{q) = ^q{j) = ^ulj. 
Explanation: If a signalling molecule diffuses from voxel j to j + 1, it means the number 
of signalling molecules in voxel j is decreased by one (hence rj{j) = — 1) and that in 
voxel j + 1 is increased by one (hence rj(j + 1) = 1). The rate at which this particular 
type of jumps takes place is proportional to the number of molecules in the j-th voxel, 
which is given by the j-th element of the state vector q. 

For convenience, we define r^ to be a zero vector and Wn{q) = 0. 

(b) For the diffusion of L from voxel j to j — 1 where 2 < j < rN-\-j{j) = — 1, rN-\-j{j — 
1) = 1 and WN-\-j{Q) = ^Q{j) = <d^L,j. For convenience, we define tat+i to be a zero 
vector and Wn-\-i{q) = 0. 

(c) For receiver 1, the vector r2Ar+i and the rate W2n-\-i{q) are used to model the forward 
reaction of the conversion of a signalling molecule L to a complex C. Specifically, 
r2N+i{Ri) = -1, r2Ar+i(A^ + l) = 1 and W2N+i{q) = ^q{Ri) = a^^l^Ri- 
Explanation: In the forward reaction, a signalling molecule is removed in the Ri-th 
voxel, hence r2N-\-i{Ri) = —1 and a complex is formed, hence r2N-\-i{N + 1) = 1 
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because the number of complexes at receiver 1 is the {N + l)-th element in the state 
vector ([s]). The rate I^2Ar+i(^) is proportional to the number of signalling molecules in 
the Ri-th voxel where the receiver is located. 

For the reverse reaction, r2Ar+2(^ + l) = — 1, ^2Ar+2(^i) = 1 and W2n-\-2{q) = k-q{N-\- 
1) = k-nc,i- 

(d) For receiver 2, the forward reaction: r2N-\-3{R2) = —1, ^2Ar+3(^+2) = 1 and W2n-\-3{q) = 
a^q{R2) = A^nL,R2' For the reverse reaction, r2Ar+4(A^ + 2) = -1, r2Ar+4(^2) = 1 
and W2N+4{q) = k-q{N + 2) = A:_nc,2- 

The RDMEX model for the 2-transmitter 2-receiver molecular communication network is 

= Y.Y.{P{q-kaATj-P{q,tmt-ta,i>) 

a=l 6=1 

J J 
+ E ^^(^ - ^o)n<l -rj.t)-Y, W,{q)P{q, t) (9) 

where 5{t) denotes the Dirac delta function. 

Let us, for the time being, assume that the first term on the right-hand side of (|9| is not there. 
If this is the case, then ([9| is of the same form as the master equation ([T]) and the equation models 
a Markov process. Given this model includes both reaction and diffusion, equation ([9| without 
the first term is known in the literature as reaction-diffusion master equation (RDME). 

The novelty of the RDMEX model is the introduction of the first term on the right-hand side 
of (|9|. This term can be viewed as a deterministic input because molecules are emitted by the 
transmitters at pre-determined times. Let us look at this term more closely. At time ^^,6, the 
a-th transmitter emits ka^h signalling molecules into the T^-th voxel (where the a-th transmitter 
is located). This means that if the system is in the state q just before the time ta^h (denoted as 
then it will be in state q + ka^h'^Ta just after ta^h (denoted as t^^). In addition, we have 
P{q^t~^) = P{q-\- ka^b'^Ta^ta b)^ which is modelled by the first term in ([9|. Note that it is possible 
to give a stochastic interpretation of see Remark [2] 

The deterministic input in (|9| can be thought of as an external arrival of the system. We will 
refer to ^ as reaction-diffusion master equation with exogenous input, or RDMEX for short. The 
name is inspired by time series models such as ARX and ARMAX [15 . 

Note that the RDMEX model is no longer a Markov process due to the deterministic arrivals. 
However, RDMEX is piecewise Markovian in the sense that, it is Markovian in between two 
consecutive deterministic arrivals. 

We have given an example of RDMEX for a 2-transmitter 2-receiver model in 1-dimension. 
The model can be readily generalised to include more transmitters and receivers. In order to 
generalise the model to 3-dimensional space, we will need to divide the space into 3-dimensional 
cubic voxels of equal volume. (The use of more complicated geometry is possible, see [16 .) The 
molecules in a voxel are only allowed to diffuse to any of its neighbouring voxels. This can also be 
readily be done. Lastly, we remark that it is also possible to use more complex receiver structure 
or to consider non-isotropic medium. 

3.2 Mean and covariance of receiver output in the RDMEX model 

We will now generalise the result of Proposition [l] to the case of RDMEX model. 

Proposition 2 For the RDMEX model in assuming that Wj{q) is a linear function of q. Let 
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^jWjiq) = Aq, then 



d{Q(t)) 



MQ{^)) + I] X] Kh^Tj{t - ta,b) (10) 



dt 

a=l 6=1 



d^{t) 
dt 



J 



= Am^mA^ ^Y.'jrjw,{{Qm (n) 



Proof: For the time evolution on the mean {Q{t)), we can start with the derivative of ([3|: 
^^^'^^^^ — ^^^^f^ and then use j9l for ^^^f^. This is fairly straightforward and uses exactly 



dt 

the same argument as the proof in [10] . 



Alternatively, one can argue the correctness of (10) as follows. Given that the difference 
between ([T]) and ([9| is the deterministic arrivals modelled by impulses, this means that between 
two consecutive deterministic arrivals, the evolution of the state in RDMEX can be described by 
a standard master equation. Hence, (|6| holds between two consecutive deterministic arrivals. It 
can be readily shown that the effect of ka^b molecules arriving at time ta,b is to add ka^b^Ta to the 



state vector. Hence the form of (10). 



For the evolution of covariance, one can follow the derivation in [W provided that the im- 
pulses are handled correctly because the multiplication of Dirac deltas (or distributions) is not 



well defined. However, one can argue the correctness of (11) using the same argument in the 
last paragraph. We know that between two consecutive deterministic arrivals, ^ holds for the 
RDMEX model. It remains to show that the covariance matrix just before a deterministic arrival 
is equal to that just after the deterministic arrival. 

Let Q{t~ij) and {Q{t~^)) be the state and mean state just before the deterministic arrival at 
time ta^b- At time t^^, just after ^^,6, the state of the system will become Q{t~^) + ka^b^r^. Also, 
the mean state at tt 



m6 



q q 

= Y.^q + KATjP{q.tlb) = mib)) + Kb^T^ (12) 

q 

Note that we have used the fact that P{q^ f^^) = P{q — ka^b^Tai'^a b) the above derivation. The 
overall result is that, at time ta,6, the state Q{t) and mean state {Q{t)) are incremented by the 
same vector ka^b^Ta- 

Given that, at any deterministic arrival, both the state and mean state change by exactly the 
same amount, therefore, deterministic arrivals do not cause discontinuity in covariance. Hence 



(11). □ 



For the rest of the paper, we will focus on studying the properties of equation (10), though we 
will present a numerical example in Section |6] to demonstrate the accuracy of (11). A detail study 



on (11) is also important and will be done in a future paper. 



Remark 2 We will now briefly discuss a generalisation of the RDMEX model. Instead of assum- 
ing a deterministic emission of exactly ka,b signalling molecules by the a-th transmitter at time 
ta^b, one may assume that the number of molecules emitted is a random variable Ka^b '^'^th mean 
(Ka^b) CL'i^d covariance cov{Ka,b)- Provided that the random variable Ka^b '^s independent of the 



state q at time ta^b or earlier, similar results to Propositioni^can be derived. For equation (10), we 
need to replace ka^b by {Ka^b), (J.nd we need to add cov{Ka^b) to the right-hand side of (11). This 
generalisation says that one can interpret ka^b ( p^ as the mean number of molecules emitted 
at time ta^b by the a-th transmitter. With this stochastic interpretation of ka^b, one can consider 
the signalling molecules are generated by an irreversible chemical reaction. 
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4 Continuum limit of RDMEX 



In section [XT] we present an example of the RDMEX model for a 2-transmitter, 2-receiver molec- 
ular communication network in an isotropic 1-dimensional medium. We also show that if the 
reaction kinetics at the receiver is linear, then the mean number of molecules in the network 



evolves according to the ODE (10). In section 4.1, we will determine the continuum limit of (10) 



as Ax 0. In order to simplify the presentation, we have so far limited our study to 1-dimensional 
but given most molecular communication networks are 3-dimensional, we will generalise the con- 
tinuum limit to 3-dimensional case as well. 

The continuum limit of the RDMEX is in fact a RDPDE. A nice property of the resulting 
RDPDE is that a closed form solution is available. This closed form solution shows that the 
output signal of a receiver can be affected by the presence of other receivers in the network. This 
will be discussed in section 14.21 



4.1 Continuum limit and generalisation to 3-dimensional space 



In this section, we will study the continuum limit of equation (10) and show that in the limit 



when the interval Ax goes to zero, (10) converges to a RDPDE and a number of chemical kinetics 
ODEs. In order that we will be able to solve the RDPDE analytically later on, we will assume 
from now onwards that the 1-dimensional medium is infinite, which means that the molecules in 
each voxel can diffuse to either of its neighbouring voxel and the state vector q is 

q = [ .... riL-2 riL-i nL,o ^l,2 .... ^c,i ^c,2 ] (13) 

where, as before, ulj is the number of L in the j-th voxel where j G Z, and nc,u is the number 



of complexes formed at the i^-th receiver. Given this state vector, we can write equation (10) as 
^^^^^^ = d((ni,,_i(i)) - 2{nL,j{t)) + {nL,,+i{t))) + ^f^^xO' - T,)ka,bS{t - ta,b) 

a=l 6=1 

2 

- ^^0' - Ru){^{nL,RAt)) - k-{nc,um Vj G Z (14) 

n=l 

^^^^^ = ^{nL,RAt))-k-{nc,uit))iovu=l,2 (15) 

where SxiJ) is the Kronecker deltaj^. 

Suppose the centre of the voxel j is at position Xj, we replace {nLj{t)) by i{xj^t)Ax where 
i{xj^t) is the mean concentration in voxel j at time t. By dividing both sides of ([l4| by Ax, 
taking the limit Ax and noting that d(Ax)^ = £), we have 



di 
dt 

d{nc,u{t)) 
dt 



dx^ 



2 

E 

a=l 



5{X - XT,a) ^ K,h5{t - ta^h) 



2 

^{x-xr^u) ± (16) 



6=1 



n=l 



dt 



k^i{xR^u,t) - k-{nc,u{t)) for u = l,2 



(17) 



where XT^a (resp. xr^u) is the centre of the voxel {Ru) where the a-th transmitter {u-th receiver) 
is located. Note that we have also used the following conversion between the Kronecker and Dirac 
deltas: limA^^o = S{x). 

This shows that in the continuum, the RDMEX converges to a RDPDE (16) and a number of 
ODEs describing the kinetics at the receivers (17). The RDPDE (16) has a simple interpretation. 



-•^Note: We use both Konecker delta and Dirac delta in this paper. They are denoted, respectively, as Sk(*) and 



S(.). 
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The second term in (16) says that the transmitter at Ta adds signalling molecules to the system 
according to time sequence ka{t)^ while the third terms says the signalling molecules are absorbed 
from the system if they form complexes at the receivers. Given that we assume that one signalling 
molecule reacts to form one complex, therefore the rate of absorption of signalling molecules is 



equal to the rate of complex formation, which is given by (17). 

Note that it is well known in literature, see [6, 16 , that a RDME with linear reaction rates 
converges to a RDPDE in continuum. In the above, we show analogues result holds for the 
RDMEX model. 

4.1.1 Generalisation to 3- dimensional space 

In order to simplify the presentation, we have so far limited to the 1-dimensional case. Given 
that most molecular communication networks are 3-dimensional, we state that in a 3-dimensional 
infinite medium the RDMEX will converge to the following RDPDE in the continuum. Here v 
denotes a point in the 3-dimensional space (i.e. v is a 3-dimensional vector) and i{v^ t) is the mean 
concentration of the signalling molecule at the location v at time t. 

2 2 

g = DSlH + Y.^{v-VT,a)K{t)-Y.5{v-VRy-^ (18) 

a=l n=l 

= k+l{xR^u. t) - k-Cu(t) for u = l,2 (19) 

where is the Laplacian in 3-dimensional space, and vr^a (resp. vr^u) is a 3-dimensional vector 
specifying the location of a-th transmitter {u-Xh receiver). Note that we have also introduced a 
new notation Cu{t) to denote the mean number of complexes {nc^u(f)) at receiver u] this is to 
simplify the notation later. 

The derivation assumes that the molecule in a voxel diffuses to a neighbouring voxel at a rate 
of d per molecule per unit time and each voxel is a cube of size A^. The parameter d is related 
to the 3-dimensional macroscopic diffusion constant I) by d = Also, the rate of formation 
of complexes at the receiver voxel is given by ^ times the number of signalling molecules L in 
the receiver voxel. Given that the derivation for the 3-dimensional is essentially the same as the 
1-dimensional case, we do not present it here. 

4.2 Solution to RDPDE 



In this section, we present a solution to the RDPDE (18) and (19), which is the continuum 
limit of the 3-dimensional RDMEX model. The key result is a closed-form expression of the 
multivariate transfer function from the transmitter signals ki(t) and k2(t) (which models the 
number of molecules injected by the transmitters into the medium at time t and can be viewed 
as the inputs to the system) to the mean number of complexes formed at the receivers ci(t) and 
C2{t) (which can be viewed as the outputs). We will divide this section into two parts. We will 
first derive the transfer function and then provide an interpretation of the transfer function. 

4.2.1 Derivation of transfer function 

The aim of this part is to derive a multivariate transfer function from ki {t) and k2 (t) to ci {t) and 



C2{t) using ([18]) and ([19]). 

We first define a few notation. Let t = a/— T, and Cu{oj)^ Ka{oj) and Liv^uo) be the temporal 
Fourier transform of, respectively, c^(t), ka{t) and i{v,t) where uj is the transform variable. It is 
shown in Appendix [A] that 

2 2 
L{V,UJ) = ^(I){V - VT,a,^)Ka{uj) -^(I){V - VR^u,^)t^^Cu{uj) (20) 
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where 



(t){v,uj) 



1 



exp(- 




(21) 



is the temporal Fourier transform of the 3-dimensional diffusion kernel - — ^^—3 exp(— Given 



(20) holds for any location we can use it to determine the concentration of the signalling 
molecules at the two receivers. By substituting v = vr^i and then v = vr^2 in ([20|, we have: 

where 



(j)ii{uj)Ki{uj) + (j)i2{uj)K2{uj) - (l)o{uj)iujCi{uj) - (j)/\,R{uj)iujC2{uj) (22) 

(j)2l{^)Ki{uj) + (j)22{(^)K2{(jO) - (l)/s^R{uj)iUjCi{uj) - (j)Q{uj)iUjC2{uj) (23) 



(/>na(^) = (/)(Vi?,n - ^T,a,^) for ii, a = 1,2 
0Ai?(^) = 0(^i?,l - ^i?,2, ^) 

= 0(0,0;) where = zero vector 



(24) 
(25) 
(26) 



One can interpret (j)ua{^) as the transfer function which models the dynamics of the diffusion 
of molecules from the a-th transmitter located at VT^a (where the molecules are injected into the 
medium) to the u-ih. receiver located at vr^w 

When signalling molecules are absorbed to form complexes, it creates a concentration gradient. 
The diffusion dynamics between the two receivers are modelled by (j)/\R{uo). Given that we assume 
that the locations of the transmitters and receivers are distinct, both (j)ua{^) and (j)/\R{uo) are well 
defined. 

The transfer function (j)Q{uj) models the local impact of absorption of signalling molecule at each 
receiver. This transfer function is unfortunately not well defined as can be seen from substituting 
V = O'm the definition of (j){v^ uo) in (21 ). The transfer function 0o(^) also appears in the modelling 
of receptor noise in chemotaxis in the biophysics literature [TT1[T2]. In fact [TTl Equation (18)] 
and pjj Equation (6)] are special cases of (18) where the input terms ka{t) are absent. Both 
[TTJ[T2] deal with the indefiniteness of 0o(^) by cutting off an integral to evaluate 0o(^) at a finite 
frequency. However, this requires us to make an assumption on the size of the receptor molecule. 
Instead, in section [sj we derive a new method to approximate (jy^iuj)^ and we will show using 
numerical examples in section |6] that our approximation gives accurate results. For the rest of this 
section, we will continue to use equations (22) and (23) with the understanding that 00 (^) is not 
well defined but can be well approximated. 



Both equations \22\ and ( |23[ ) are obtained from ([18|. We still need to work on ( |19| ). By taking 
the Fourier transform of (19), we have 



Pi(^) 

P2(^) 



L{vR^i,uj) 

C2{UJ) 
L{vr^2,^) 



(27) 
(28) 



where pi{uj) and P2{^) are transfer functions that model the reaction kinetics at the receivers. 
Given that we have assumed that the binding and unbinding rate constants at both receivers are 
identical, it is not surprising that pi{oj) and p2{^) are the same here. It is straightforward to 
generalise to the case where the receivers have different reaction kinetics. 

By using equations (22), (23), (pTl and (28), we can eliminate L{vr^i^uj) 



and L{vr^2-,^) to 



obtain the transfer function from the inputs ki{t) and k2{t) to the outputs Ci{t) and C2{t)\ 



( 



C2 



" pi 







00 0Ar 





P2 _ 




0Ar 00 



Pi 







P2 



hi 



n2 

/>22 



K2 



(29) 
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where, for conciseness, we have dropped the dependence on transform variable uo. Equation (29) 



is the key result of this section. It is the solution to the RDPDE (18) and (19), which are in turn 
the continuum limit of the mean concentration in the stochastic RDMEX model. We will provide 
some physical interpretation of (|29l 



in a moment. Before that, we want to point out that (29) 
can be used to compute Ci (t) and C2 (t) given ki {t) , k2 (t) and the system parameters by numerical 
Laplace transform. This will be done in Section [6j 



Note that it is numerically more efficient to solve for the receiver outputs using (29) rather 



than (10). This is because one also needs to solve for the number of signalling molecules in the 



voxels in (10) but this is not needed when (29) is used. Since the number of voxels is much larger 



than the number of transmitters and receivers, numerical solution via (29) is more efficient. 



4.2.2 Interpretation of transfer function (29) 



Equation (29) may not look easy to interpret in the first instance, so we will first specialise it to 
the case of 1-transmitter and 1-receiver. In this case, we have 

Pi^ii 



1 + iUJ(t)Qpl 



(30) 



One can readily show that the input-output transfer function in ( 30 ) corresponds to the block 
diagram representation of figure |3] The negative feedback loop occurs because the net number 
of signalling molecules at the receiver i{vR^i^t) is given by the difference between those that 
arrive via diffusion (modelled by the feedforward block of (pn) minus those reacted to form the 
complexes (modelled by the feedback block of LUJ(j)o). The forward loop consists of which 
models the diffusion dynamics of signalling molecule from the location of the transmitter to that 
of the receiver, and pi which models the conversion of the signalling molecules to complexes. 

We will now take a closer look at the denominator of (30) with the aim to determine the 
strength of the feedback term iuj(j)Qpi. It has been shown in [111, 1X2] that 0o oc at low 

frequency. Consider the case where D ^ kj^ or ^ ^ which corresponds to the situation where 
the chemical kinetics is not diffusion-limited. We have LUJ(j)opi ~ base on the expression of pi 
in (27) and consequently Ci ~ picj^iiKi. Since the chemical kinetics is not diffusion-limited, the 
chemical kinetics and diffusion are basically " decoupled" , so the transfer function from Ki to Ci 
is the multiplication of the transfer function from Ki to L{vr^,uj) (which models diffusion) 
and the transfer function pi from L{vr^^uo) to Ci (which models reaction kinetics). We will show 
in section 6 using numerical examples to show that the transfer function ^ Pi^^ii holds 

>k+. 



when D ^ 

An end-to-end frequency model for molecular communication with one transmitter and one 
receiver has earlier been derived in [13]. The model assumes decoupling of diffusion and chemical 
kinetics [13, Equation 1] which means that the model in flT holds when diffusion occurs at a 
faster rate than reactions. 

Given the interpretation of the 1-transmitter 1-receiver case as a system with feedback in figure 
[sj it can be shown that equation ( 29 ) corresponds to multivariate feedback system with 2 inputs 
and 2 outputs. The block structure of the multivariate feedback system is the same as that in 
figure [3] but we need to replace the single-input single-output transfer functions ^n, pi and 0o by 
their multivariate counterparts of ^i, 7^ and <l>o in (29). 



For the 2-transmitter 2-receiver case, we see from (29) that the response at each receiver is 
affected by both transmitters, as well as by the other receiver. Let us assume that the transmitters 
and receivers form two communication pairs where transmitters 1 and 2 intend to communicate 
with, respectively, receivers 1 and 2. We want to determine the unintended signal at the receivers. 
Let us assume for the time being that the two receivers are sufficiently far apart so that (j)Ar{^) 
is negligible compared with ^o(^)- Iii addition, we assume the two receivers are identical, so 
pi{oj) = p2{^)- In this case, we can simplify (29) to 



Ci 

C2 



Pi 



1 + LCOCpoPi 









4'2l (t>22 




. ^2 . 



(31) 
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Comparing with the single-transmitter single-receiver transfer function in (29), we can see that 
the unintended signal due to transmitter 2 on receiver 1 is i^^^^^2- One can readily see that 
the magnitude of this unintended signal can be reduced if transmitter 2 is well separated in space 
from receiver 1 because (/) is a decreasing function of distance. Similar conclusion can be drawn 
for transmitter 1 and receiver 2. Note that the above argument requires that the receivers are well 



separated. If this is not the case, the matrix inverse in (29) will create a complicated mixture of 
signal at both receivers. 

In general, spatial separation is a strategy to reduce the magnitude of unintended signal that 



one communication pair has on the others. It is interesting to point out that (29) allows one to 



explore other methods to reduce the magnitude of the unintended signal. An interesting case to 
study is if the transmitters emit molecules at different frequencies and the receivers are frequency 
sensitive. We will not explore this further here and leave this for future work. 



5 Discrete solution for mean concentration in RDMEX 



In the last section, we presented a closed- form solution to the RDPDE ([Ts]). A problem that 



we face is that the frequency response ^o(^) in (26) is not well defined. This problem arises 
because the size of voxel, in the continuum limit, becomes zero. Therefore, a solution to overcome 
this problem is to consider finite voxel size instead. This means that we need to work with the 



3-dimensional analogue of equations (14) and (15). 



Consider an isotropic 3-dimensional space. We divide the space into identical cubic voxels of 
volume each. We index the voxel using a 3-dimensional vector ^ = [i j k] where i^j^k G Z; 
note that we will use ^ and [i j k] interchangeably in this section. We use ^T,a (^i?,n) to index the 
voxel that the a-th transmitter {u-th receiver) is located. We will also use intuit) to denote the 
mean concentration of signalling molecules at the voxel where the u-th receiver is located. Let 
Cu{t) denote the mean number of complexes at the u-th receiver for this discrete model and Cu{^) 
be its continuous Fourier transform. (Just to avoid any possible confusion. We discretise only 
space, not time. So, t remains continuous.) 

Let iij^k{t) denotes the mean concentration of the signalling molecule in voxel i^j^k. The 
mean concentration in a voxel is given by the mean number of molecules divided by the volume 



of a voxel which is A^. One can show that the generalisation of (14) and (15) to 3-dimensional 
space — with mean concentration of signalling molecules, rather than mean number, per voxel 
is: 



dt 



dcujt) 
dt 



a=l u=l 

fc+€_R,„(t) - k-Cu{t) for u = l,2 



dt 



(32) 
(33) 



where D = 



One can readily show that the continuum limit of i32h and (33) is (18) and il9\ 



Let Ldi^^uj) denote the temporal Fourier transform of iij^kit)- (Recall that ^ = [z, j, k].) It is 
shown in Appendix [B] that 



E ~ ^T,a, ^)^a(^) - E ~ ^)^^^n(^) 



(34) 
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where 



1 



\k\+l 



(35) 



where Wx and Wy are complex variables and the contour C is the unit circle on the complex plane; 
also, Wz^ is the solution of the following quadratic equation in Wz with modulus less than unity: 



(2 + (W, - W-'Y + (Wy - Wy-'Y + uo — )W, + 1 = 



(36) 



Note that (34) is the discrete space analogue of (20). Thus one can identify (j){0^uj) with 
tlj{0^uj). In addition, one can show that: 



C2 



= I 



iUO 



Pi 











r 


pi 









i^i2 









P2 _ 






^0 







P2 . 




_ i^2i 


^22 




. ^2 _ 



(37) 

^T,a,^), 



where for conciseness we have dropped the dependence on cj, and il^uai^ 
V^Ar(^) = i^{^R,i -^i?,2,^) and V^o(^) = V^(0,a;). 

Numerical integration can be used to compute ?/;(^,a;) in (35). This will be used in Section [6j 



6 Numerical examples 
6.1 Overview 



In this section, we will give a number of numerical examples to show that equations (29) in section 
4.2 and (37) in section [s] can be used to accurately predict the mean output of the receivers 
in the stochastic model RDMEX. We will also use these numerical examples to illustrate the 
issues of using molecular signalling for communication. We will present three sets of results: 



single-transmitter single-receiver in section [6^ single-transmitter two-receiver in section [6^ and 
two-transmitter two-receiver in[ 



In order to verify the accuracy of (29) and (37), we will use simulation to compute the mean 



receiver output. One method is to simulate RDMEX many times and compute the mean. Alter- 
natively, one can use the fact that if the number of molecules is large, then the behaviour of one 
simulation run is fairly close to the mean. We will mainly use the latter method in this paper. 
We simulate the RDMEX model using the r-leaping method [17J. The r-leaping algorithm uses a 
constant time step to advance the simulation and is a faster alternative to the Gillespie's algorithm 
[17]. We will refer to the simulation result as RDMEX-r. 

The fluidic medium is assume to have a D of 0.05. (Since the parameters in the diffusion- 
reaction system can be scaled to some dimensionless quantities \1S^ Section 8.2], we do not specify 
the units for the parameters here.) The locations of the transmitters and receivers, as well as the 
reaction rate constants, vary between experiments and will be specified later. Different transmitter 
signals will be used to demonstrate the accuracy of the RDMEX model. Our goal is to compare 
the output ci{t) and C2{t) from RDMEX-r with that from the following analytical models: 



1. Equation (|37|). We wih refer to this as RDMEX-M. 
2 



Equation (29) of the continuum model with (/)(0,co') replaced by ^/^(O,^^). We will refer to 
this as RDPDE. 



3. Equation (29) of the continuum model with (/)(0,co') approximated by ^^^/^ • This approxi- 
mation is inspired by the one used in [Til Il2i where we have replaced the size of receptors 
by the voxel dimension parameter A. We will refer to this as RDPDE-X. 

4. For the single-transmitter and single-receiver case, we consider the "decoupled" model Ci = 
Pi^ii^i- This will be referred to as DE. 
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For all these analytical models, we first determine the Laplace transforms Ci and C2 by us- 
ing the Laplace transforms Ki and and the transfer function. We then invert the Laplace 
transform numerically using the matlab function invlap.m [19^ . 

6.2 Single-transmitter single-receiver case 

The system consists of a transmitter at the voxel [0,0,0] and a receiver at voxel [3,0,0]. 



6.2.1 Model accuracy and the effect of /c+ 

In this set of experiments, the transmitter emits 10 molecules every 10~^ time units for a duration 
of 0.2 time units and then stops emitting for 0.3 time units. The signal ki{t) is obtained by 
concatenating this emission pattern 3 times. The value of k- is 0.05. We determine the output 
signal ci(t) for t e [0,2]. 

We use two different values for Figure |4] shows the results for /c+ = 2.5 x 10~^. The figure 
compares the mean number of complexes formed in the time interval [0,2]. Results are obtained 
from RDMEX-r (simulation) and RDPDE, RDMEX-M, RDPDE-X and DE (analytical models). 
It can be seen that both RDPDE and RDMEX-M (our analytical solutions) match RDMEX-r 
well. The model RDPDE-X does not give good approximation because the voxel size is not a good 
approximation for the receptor size. Henceforth, we will not consider RDPDE-X further. The 
decoupled model DE does not match RDMEX-r for this value of 

We then change to 2.5 x 10""^. The results are plotted in Figure [s) Both RDMEX-M and 
RDPDE again match RDMEX-r well. We also see that the decoupled model DE gives better 



prediction than before. This validates the discussion in Section |4.2.2| that the decoupled model 
holds when /c+ is sufficiently small. 



6.2.2 Accuracy of mean and standard deviation computation 



In this experiment, we use a different transmitter signal to show the accuracy of RDMEX-M and 
RDPDE. We define two transmitter symbols sq and Si. The symbol duration is 2 time units. 
When a transmitter sends 5i, it emits 10 molecules every 10~^ time units for a duration of 0.2 
time units and then stops emitting for 1.8 time units. When a transmitter sends sq, it does not 
emit any molecules for 2 time units. The transmitter signal kiit) in this experiment is simply si. 
The receiver parameters are /c+ = 2.5 x 10~^ and k- = 8. We determine the output signal ci(t) 
for t e [0,2]. 

We first verify the accuracy of using RDMEX-M and RDPDE to compute the mean of receiver 
output. For this experiment, we simulate the system using RDMEX-r 125 times and compute the 
mean receiver output as the reference. Figure [6] shows the mean receiver output from RDMEX-M, 
RDPDE and RDMEX-r. It can be readily seen that RDMEX-M and RDPDE are accurate also 
for a different transmitter signal. 



Our next goal is to verify the accuracy of using equation (11) to determine the standard 
deviation of the receiver output. We extend the RDMEX-M to solve for both the receiver output 
as well as the system states, which are the number of signalling molecules in the voxels. The 



system states are used as the input to (11) and numerical integration is used to solve for (11). For 



verification, we simulate the system using RDMEX-r 125 times to compute the standard deviation 
of the receiver output. Figure 7[plots the standard deviation of the receiver output from the two 
methods. It can be seen that (H]) is accurate. 

One can envisage using the symbols si and sq to encode the communication between the 
transmitter and the receiver. The communication scheme is similar to ON-OFF keying. For 
decoding, the receiver can use, say, the peak number of complexes to detect the symbol transmitted. 
If the peak number of complexes is above a threshold, then si has been transmitted; otherwise, 
So has been transmitted. 



Remark 3 The above method of calculating the variance using (11) is computational intensive. 



A future work is to derive efficient algorithm to solve (11). 
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6.2.3 Assumption of linear receiver model 

In this paper, we assume that the rate of formation of complex C is a hnear function of the number 
of signahing molecules in the receiver voxel, according to the following chemical reaction kinetics: 



ki- 



(38) 



Under certain assumptions, the above reaction kinetics can be used to approximate more 
complex reactions. Consider the following chemical reactions: 



L + E^I 

gi- 

C ^ L 



C + E 



(39) 
(40) 



type where a molecule L reacts with an enzyme E to form 
following by the decomposition / into a product C and the enzyme 



Reaction ( 39 ) is of Michaelis-Menten 
an intermediate product / 
E. It can be shown that, for suitable choice of reaction constants in (|39|), the reaction kinetics of 



(39) can be approximated by the forward reaction in (38). Similarly, reaction (40) can be made 



identical to the reverse reaction in ( 38 ) by choosing g^, 



We want to show that a receiver with reactions (39) and (40) gives similar output signal 



compared to one with reactions (38). We implement the receiver kinetics (39) and (40) in r-leaping 



simulation and refer to this method as RDMEX-r-MM. We compare this against RDMEX-M with 
identical parameters to those used in section 6.2.2 The parameters in reactions ([39l) and (40) 



have been chosen to approximate those in (38). The results are plotted in Figure [s] and it can be 
seen that RDMEX-M is able to approximate more complicated reactions. 



6.3 Effect of A 

In this section, we study the effect of the size of voxels. We use the same parameter setting for 
the single-transmitter single-receiver case in Section 6.2.1 We assume both the transmitter and 
receiver is a cube whose length of an edge is x = 0.01. We use four different voxel sizes with 
^ = f 7 3 arid ^. Since the size of the transmitter (or receiver) is a constant, this means that 
the transmitter (or receiver) occupies, respectively, 1, 8, 27 and 64 voxels for these 4 different voxel 
sizes. We assume that the emissions from the transmitter is uniformly distributed across all the 
voxels that it occupies. Figure |9] show the mean number of complexes at the receiver, which is the 
sum of the number of complexes at all receiver voxels. It can be seem that the predicted output 
for A = |, I and ^ are almost the same. This shows that as long as A is sufficiently small, the 
prediction is independent of A. 



6.4 Single-transmitter two-receiver case 



Equation ( [29| ) shows that when there are multiple receivers, it is possible for a receiver to affect 
the output of another receiver. We will illustrate this phenomenon. We consider three different 
networks. Network consists of a transmitter and two receivers, called 1 and 2. The transmitter, 
receivers 1 and 2 are located, respectively, at voxel [0,0,0], [1,0,0] and [2,0,0]. Network 1 is 
composed of the transmitter and receiver 1 of network 0. Network 2 consists of the transmitter 
and receiver 2 of network 0. For all the three networks, the transmitted signal is 5i and the 
receiver parameters are /c+ = 2.5 x 10~^ and k- = 8. 



Figure [To] shows the output for receiver 1 for networks and 1, while Figure pT] shows the output 
for receiver 2 for networks and 2. (The curve with labelled ltl5r in Figure flT| will be explained 
later.) It can be seen that the output of receiver 1 is almost the same for both network (receiver 
2 present) and network 1 (receiver 2 absent). However the output of receiver 2 for network 
(receiver 1 present) is very different from that in network 2 (receiver 1 absent). Specifically, when 
receiver 1 is present, the output of receiver 2 has a lower peak number of complexes and a higher 
number of complexes at the tail. 



15 



An explanation of how receiver 1 affects the output of receiver 2 is as fohows. Note that 
receiver 1 is situated in between the transmitter and receiver 2. Some signahing molecules that 
reach receiver 2 have to pass through the voxel containing receiver 1. When these signalling 
molecules are in the receiver 1 voxel, some of them react to form complexes and are held up in 
the voxel. This means less signalling molecules reach receiver 2 in the early part of the symbol 
duration, thus resulting in a lower peak number of complexes in receiver 2. During the later part 
of the symbol duration, the complexes in receiver 1 dissociate to release the signalling molecules. 
Some of these signalling molecules, which are held up earlier in receiver 1, reach receiver 2 later 
on. This means more signalling molecules reach receiver 2 in the later part of the symbol duration 
if receiver 1 is present. This explains the behaviour at the tail of the output of receiver 2. 

When the number of receivers in a network is high, the output of some receivers in a network 
can be highly affected by the presence of the other receivers. We create a network with a large 
number of receivers by adding an additional 14 receivers to network 2, making a total of 15 receivers 
in the network. We again focus on the output of receiver 2. The curve labelled with ltl5r in 



Figure 11 shows the output of receiver 2 for this 1-transmitter 15-receiver network. It can be seen 
that the receiver output in this network is very different from that receiver 2 in network 2. The 
reason is that the other 14 receivers are affecting the output of receiver 2. 

The above results mean that the design of molecular communication need to take all receivers 
in the network into consideration. For example, if receiver 2 uses a threshold on peak number of 



complexes to detect si, we can see from Figure 11 that a threshold that works for the 1-transmitter 



2-receiver network may not necessary work for the 1-transmitter 15-receiver network. 
6.5 Two-transmitter Two- receiver case 

We consider a molecular communication network with 2 transmitters located at [0, 0, 0] (trans- 
mitter 1) and [4,0,0] (transmitter 2), and 2 receivers located at [1,1,0] (receiver 1) and [2,0,1] 
(receiver 2). Both transmitters use si and sq defined earlier as their transmission symbols. The 
signals for transmitters 1 and 2 are, respectively, sqSi and siSq. 



We first verify the accuracy of our proposed models. Figures [12] and [13] show the output signals 
for, respectively, receivers 1 and 2. The three curves in each figure are obtained from RDMEX-r 
(simulation) and RDPDE and RDMEX-M (our analytical solutions). It can be seen from both 
figures that the prediction from both RDPDE and RDMEX-M match that of RDMEX-r. In fact, 
the curves in the figures match so well that they overlap. 

We now assume the transmitters and receivers form two unicast communication pairs: trans- 
mitter 1 communicates with receiver 1 while transmitter 2 communicates with receiver 2. If 
transmitter 1 were the only transmitter, then receiver 1 should have a zero signal in the first 
symbol duration (between and 2 time units). However, Figure 
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shows that receiver 1 has a 
non-zero signal during the first symbol duration due to the transmitter 2 sending an si during 
this time. During the second symbol duration, the output signal of receiver 1 is due entirely to 
transmitter 1. If receiver 1 uses a threshold based detector, then a suitable choice of threshold 
will enable receiver 1 to correctly decode the two symbols sent by transmitter 1. 



Let us now consider the output signal of receiver 2 shown in Figure [13] The signal in the first 
symbol duration is due to transmitter 2 (the intended signal) while that in the second symbol is 
due to transmitter 1 (the unintended signal). If receiver 2 uses a threshold based detector, then due 
to the unintended signal, a bit error will occur in the second symbol duration. This is a typical 
example of bit error when multiple transmitters and receivers share a common communication 
channel. 



7 Related work 

Molecular communication networks can be divided into two categories, according to whether 
they are natural or synthetic. Natural molecular communication networks are prevalent in living 
organisms. Their synthetic counterparts, though still rare, do exist. For example, [21 presents 
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a system with multiple genetically engineered cells that use cell signalling to coordinate their 
behaviour. 

The modelling of natural and synthetic molecular communication networks is studied in dif- 
ferent disciplines. The former is mainly studied in biophysics and mathematical physiology, while 
the latter in synthetic biology. There is also a recent interest in the engineering community to 
study molecular communication networks from a communication theory point of view [13 . 22 , 23^. 
This gives rise to a new research area called nano communication networks [1 . 

Despite the fact that molecular communication networks are studied in diverse disciplines, the 
set of mathematical models that are being used are similar. This is not surprising given that the 
primary goal is to model diffusion and reaction kinetics. The classes of mathematical models being 
used include molecular dynamics, master equation, partial differential equation (PDE), Fokker- 
Planck equation, Langevin equation and others [6]. We will focus on the first three classes of 
models in this discussion. 

Molecular dynamics is commonly used in simulation of molecular communication networks. 
Many examples of simulators exist, especially for natural molecular communication networks, see 
[24] for a recent overview. For synthetic networks, a recent example is |25]. By analysing the 
molecular dynamics of transmitters and receivers, [26 characterises the noise in transmitters and 
receivers as, respectively, sampling and counting noise. 

There are ample examples in using PDE — in particular diffusion PDE, telegraph equation 
and RDPDE — to model molecular communication. For natural networks, [TTJ [12] use RDPDE 
to study the noise in receptor binding in chemotaxis, and [^ uses RDPDE to study signalling 
cascades. However, these papers do not consider the transmitters. For synthetic networks, tele- 
graph or diffusion PDEs (or their kernels) have been used to characterise the diffusion of signalling 
molecules in ^| E] EH [22l [23] and others. However, these papers do not consider the coupling 
effect between diffusion and receiver reaction kinetics. In our earlier work in [30 , we use a RDPDE 
in the form of (18), as a deterministic model for molecular communication network. The RDPDE 
in [30 is solved numerically and no analytic solution is provided. In this paper, we derived a RD- 
PDE model (18) for molecular communication and provide an interpretation of the model as the 
mean receiver output of molecular communication networks. In addition, we present an analytical 
solution to this RDPDE and show that it can be used to accurately predict mean receiver output 
in molecular communication networks. 

For some time, RDME has been considered to be a phenomenological model because it diverges 
in certain cases [31]. Fortunately, the problem has been resolved in [32 and there is now a firm 
theoretic basis for RDME. There are many examples of work that use RDME to model natural 
molecular communication networks, see [TO, "33]. However, these papers do not consider the 
transmitters. The use of RDME in studying synthetic molecular communication networks appear 
to be novel. To the best of our knowledge, our RDMEX model, which is formed by coupling time 
sequences of signalling molecule emission pattern with RDME, has not been proposed before. The 
proposed RDMEX model is one of the novel contributions of this paper. 



8 Conclusions and future work 

In this paper, we have proposed a new stochastic model called reaction-diffusion master equation 
with exogenous input (RDMEX) for modelling molecular communication networks with multiple 
transmitters and multiple receivers. We show that we can readily derive the mean and covariance 
of receiver output of RDMEX model for the case where reaction kinetics at the receiver is linear. 
We then solve the mean receiver output of RDMEX model using two different methods. In the first 
method, we derive the continuum limit of RDMEX and present a closed-form expression of the 
solution to the resulting reaction-diffusion partial differential equation. In the second method, we 
solve the mean receiver output of RDMEX explicitly by using transform techniques. We present 
numerical examples comparing the accuracy of our analytical solutions against simulation. We 
find our analytical solutions give accurate prediction of mean receiver output of molecular com- 
munication networks. This paper has focused on studying the mean receiver output of molecular 
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communication networks. However, one needs to understand the properties of noise in these net- 
works in order to evaluate their performance, such as bit error rate or capacity. Future work 
includes the study of the properties of noise in the RDMEX model and using the RDMEX model 
to evaluate the performance of molecular communication networks. 



A Derivation of ([20|) 



Equation ( |2Q| ) can be derived by using a couple of different methods, e.g. Fourier transform or 
Green's function. We will use Green's function here. 



Note that the RDPDE (18) is an inhomogeneous partial differential equation where the last 
two terms act as the forcing function. Let G{v^ t) denote the three dimensional kernel of the linear 
diffusion equation || = DV^l. Let also ^s,t denote convolution in space and time, and *t denote 
convolution in time only. By using the theory of Green's function, we have 



dcu{t) 



£{v,t) = G{v,t)^s,T<Y,S{v-VT,a)ka{t)-Y,S{v-VR,u) 
\.a—l u=l 
2 2 
= ^ G{V - VT,a, t) *T K{t) - ^ G{V - VR^u.t) *t 



dt 



(41) 



Now we can obtain (20) from applying temporal Fourier transform to (41), noting that (j){v,uj) 
is the temporal Fourier transform of G{v^t). 



B Derivation of (34) 



Note that (p2|) is a linear difference equation in space and a continuous differential equation in 



time. We can draw a parallel with the theory of Green's function and note that (34) holds if 
ip{£^^uj) is the temporal Fourier transform of the solution of the following equation: 



dt 



(42) 



This equation can be solved by using Fourier transform on the continuous variable t and z- 
transform on the discrete variables j and k. Let Hi^j^k(uj) denote the Fourier transform of 
hij^k{t)' The multi-dimensional z-transform of Hij^ki'^) is defined as: 



oo oo 



i= — ooj= — oc k= — oo 



By using (42), it can be shown that 



^iW,,,Wy,W,,LO) 



DA{{W^ - W^^Y + {Wy - Wy^Y + {W, - w^^Y - 



We first compute inverse z-transform with respect to Wz- We re- write (44) as 

-1 W, 



DA {Wz-Wz,){Wz-Wr*') 



(43) 



(44) 



(45) 



where VF^*, as defined in section [sj is the root of ( [36| ) with modulus less than unity. Note we have 
also used the fact that the roots of (|36]) are reciprocal of each other. It can now be shown that 
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the inverse z-transform of "^{Wx^ Wy^ Wz^uj) with respect to Wz is 



-1 wi'J^' 



(46) 



DA Wl - 1 

If we apply the standard inverse z-transform contour integrals, with respect to Wx and Wy^ to 
this expression, then we obtain the formula for ip{v^uj) given in (35). 
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Figure 1: An illustration of a molecular communication network with two transmitters and two 
receivers. The signalling molecules diffuse in a fluidic medium and may bind with the receptors 
at the receviers. 
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Figure 2: The voxels in the 1-dimensional RDMEX model. 
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Figure 3: Block diagram showing the transfer function from the input ki{t) to the output ci{t) of 
a molecular communication network with one transmitter and one receiver. 
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Figure 4: The mean number of complexes in the 1-transmitter 1-receiver network for K = 2.5 x 
10-^ (Section 
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Figure 5: Th e mean number of complexes in the 1-transmitter 1-receiver network for K = 2.5 x 
IQ-^. (Section |6.2.1[) 
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Figure 6: The mean number of complexes in the 1-transmitter 1-receiver network with transmitter 



symbol si. (Section 6.2.2) 
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Figure 7: The standard deviation of the number of complexes in the 1-transmitter 1-receiver 



network with transmitter symbol Si. (Section 6.2.2) 
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Figure 8: The mean number of complexes in the 1-transmitter 1-receiver network for section [6. 2. 3[ 
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Figure 9: This figure shows the effect of the size of a voxel on the mean number of complexes. 
The A values used are f ? f f- 
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Figure 10: This figure compares the output signal of two different networks. Network consists 
of a transmitter and 2 receivers (receivers 1 and 2). Network 1 consists of the transmitter and 
receiver 1 of Network 0. The figure shows the output signal of receiver 1 for these two networks. 




Figure 11: This figure compares the output signal of three different networks. Network consists 
of a transmitter and 2 receivers (receivers 1 and 2). Network 2 consists of the transmitter and 
receiver 2 of Network 0. The third network (label as ltl5r) consists of network 2 plus another 14 
receivers. The figure shows the output signal of receiver 2 for these networks. 
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Figure 12: The output signal of receiver 1 in the 2-transmitter 2-receiver network. 




Figure 13: The output signal of receiver 2 in the 2-transmitter 2-receiver network. 
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